NPSS9  -34-001 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


USE  OF  THE  TENSOR  PRODUCT  FOR  NUMERICAL  WEATHER 

PREDICTION  BY  TEE  FHTITE  ELEMENT  METHOD  -  PART  1 

R.  E.  NEWTON 
April' 1984 


Final  Report  for  Period 
October  1932  -  September  1933* 


Approved  for  Public  Release;  Distribution  Unlimited 
>ared  for: 

D  208.14/2     ■!  Environmental  Prediction  Research  Facility 
NPS-69-84-001  ere^'  California  93943 


' 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


CoomDdore  R.  H.  Shumaker  ?•  A-  Schrady 

Superintendent  ***** 


The  work  reported  herein  was  supported  by  the  Naval  Environmental  Predic- 
tion Research  Facility. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 
This  report  was  prepared  by: 


unclassified 


SECURITY   CL* 


.SlFlCATIOM   OF   THIS  PAGE  rt»7n>n  Dtilm  Rntarad) 


REPORT  DOCUMENTATION  PAGE 


t.     REPOHT   NUK3£H 

NPS69-34-001 


£Y  KNOX  LIBRARY 

GRADUATE  SCHOOL 
• 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  KORM 


2.   GOVT    ACCESSION   NO 


«.     TITUE  (and  Subtltla) 

USE  OF  THE  TENSOR  PRODUCT  FOR  NUMERICAL  WEA3EER 
PREDICTION  3Y  THE  FINITE  ELEMENT  MEEKD  -  PART  1. 


7.     AUTHORf*; 

R.  E.  Newton 


3       RECIPIENT'S  CATALOG   NJWflEH 


S.     TYPE   OF    REPORT   h    PERIOO   COVERED 

Final,  October  32  - 
September  83 


ft.  PERFORMING  ORG.  REPORT  NUMScS 


i.     PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California  93943 


S.     CONTRACT   OR   GRANT   NUM8ERC*,) 

N668568AWR84103 


1C.     PROGRAM   ELtMCNT.  PROJECT.   TASK 
ANEA   4   WOHK   UNIT   NUMUEMS 


II.     CONTROLLING  OFFICE  NAME   AND  ADORESS 

Naval  Environmental  Prediction  Research  Facility 


Monterey,  California  93943 


14.     MONITORING  AGENCY  NAME  &    ADORESSf//  diltarant  Irom  Controlling  Oltlca) 


611531-1 


12.     REPORT    DATE 

April   1984 


13.     NUMBER  OF  PAGES 


11.     SECURITY   CLASS,  foi  rhla  rorort) 

Unclassified 


\la      DP-CLASSIFICATION    DGWNGRAOi.-lG 
SCHEDULE 


16.     DISTRIBUTION  STATEMENT  (ot  thta  Hapott) 

Approved  for  Public  Release;  Distribution  Unlimited 


17.     DISTRIBUTION  STATEMENT  (of  tha  ab.tr.ei  antarad  In  Block  10.  H  dltlarant  from  Report) 


18.     SUPPLEMENTARY  NOTES 


"»""   KEY  WORDS  (Conttnua  on  rarataa  ml  da  II  n.c...«y  and  tdantlly  by  block  numbar) 

Finite  element,  numerical  weather  prediction,  tensor  product 


20.     ABSTRACT  (Conttnua  on  tavataa  alda  H  naca*o**y  and  Idamity  by  block  numbar) 


This  report  concerns  application  of  the  tensor  product  in  solving 
large  sets  of  linear  equations  arising  in  Numerical  Weather  Pre- 
diction problems.   The  coefficient  matrix  considered  is  called 
the  "mass"  matrix  in  finite  element  parlance.   Theory  of  the 
tensor  product  resolution  is  presented,  along  with  methods  for 
imposing  Dirichlet  and  cyclic  boundary  conditions.   Operation 


/ 


DO     IM,nS     1473  roiTIOM  OF    1  MOV69  ISOKSOLCTe 

S/N    0102-014-  6*01   I 


unclassified 


SKCUmTY   CLASSIFICATION   OF    .MIS   PAGC  ("fan  Data   in,a(ad) 


unclassified 


(utuWTv  Cl  A»«|g'C*  tiqm  Q*  This  »«fltfWmi  n«u  «»f»»rf- 


20 


A3 S TRACT 


counts  and  storage  requirements  are  compared  with  corresponding 
numbers  for  some  widely-used  solution  algorithms.   Listings  of 
Fortran  programs  for  implementing  the  tensor  product  solution 
system  (TENSOR)  and  testing  it  are  given. 


1473 


DD     Form 
1  Jan  «3 
S/N    0102-014-6601 


unclassified 


$tCU<»lTV  CLAMI'ICATIOM  OW  TMI»  »*OC/"»«"  ©•*•  ««'•'•*» 


lV™°*  UBRARV 


USE  OF  THE  TENSOR  PRODUCT  FOR  NUMERICAL  WEATHER 
PREDICTION  BY  THE  FINITE  ELEMENT  METHOD  -  PART  1. 

Introduction 

In  Ref.  1  Hinsman  has  developed  a  Finite  Element  (FE) 
program  for  Numerical  Weather  Prediction  applications.  The 
grid  employed  is  rectangular  with  nodes  at  the  intersections 
of  north-south  and  east-west  lines.  It  was  shown  by  Stani- 
forth  and  Mitchell  (Ref.  2)  that  the  coefficient  matrices 
for  such  a  grid  could  be  expressed  as  tensor  products.  In 
these  products  the  factors  are  matrices  which  depend  solely 
on  grid  spacing  in  the  two  orthogonal  directions.  This 
report  deals  with  the  coefficient  matrix  called  the  "mass" 
matrix  in  FE  parlance.  (In  Refs.  3  and  4  applications  of 
the  tensor  product  resolution  to  the  FE  "stiffness"  matrix 
are  considered. )  The  theory  which  underlies  the  economical 
computational  scheme  based  on  the  mass  matrix  resolution  is 
first  presented.  Next,  the  number  of  floating  point  opera- 
tions and  the  number  of  storage  locations  needed  for  the 
coefficient  matrix  of  this  scheme  are  compared  with  those 
required  by  other  better-known  algorithms.  A  set  of  FORTRAN 
subroutines  for  implementing  the  tensor  product  scheme 
(TENSOR)  is  given  in  Appendix  B. 

Theory 

Consider  the  grid  shown  in  Fig.  1.   There  are  n  spaces 
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n  columns 
Fig.    1.      Node  numbering   and   spacing. 

along   each  of   e   grid    lines    in  the   east-west    direction.      Node 
numbering   is    from  west    to      east    along   successive   grid    lines, 
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beginning  in  the  southwest  corner.  There  is  a  cyclic  bound- 
ary condition  in  the  east-west  direction  so  that  the  node 
number  appearing  at  the  beginning  of  each  horizontal  row  is 
repeated  at  the  end.  Spacings  of  the  horizontal  and  the 
vertical  grid  lines  are  not  necessarily  uniform. 

The  computational  problem  addressed  here  is  the  solution 
of  the  equation 

Mw  =  v  <1> 

where  M  (the  "mass"  matrix)  is  a  square,  symmetric  matrix  of 
size  ne  and  w  and  v  are  column  vectors  of  height  ne.  M  and 
v  are  input  quantities  and  w  is  sought.  The  tensor  product 
representation  of  M  is 

M  =  MB  *  MA  <2> 

where  MB  is  a  square,  symmetric,  tridiagonal  matrix  of  size 
e  and  MA  is  also  square,  symmetric  and  of  size  n.  MA  is 
tridiagonal  except  for  nonzero  elements  in  upper  right  and 
lower  left  corners.  MB  depends  solely  on  the  north- south 
node  spacing  b  and  MA  depends  upon  the  east -west  node  spac- 
ing a  .  The  asterisk  (*)  denotes  the  tensor  product. 
Explicit  expressions  for  matrices  MA,  MB,  and  the  tensor 
product  are  given  in  Appendix  A. 

Let  MB  be  represented  as  (e  =  3) 

mbn     mbi2       0 
MB     =      mb2i     mb22     mb23 
0        mb  3  2     mb  3  3_ 
If  we  partition  w  and  v   into   e   n  x   1   subvectors    so    that 


w  = 


wi 
wn 

V    = 

h  1 

_wm_ 

_vi". 

<4> 


we  may  use   <2>,    <3>   and   <4>   to   rewrite   <1>    as 
mbn  MA  w.  +  mbi2  MA  w,j  =  Vj 

mb21  MA  w.  +  mb22  MA  Wjj  +  mb23  MA  w^j  =  vn 


mb 


32  MA  w.j  +  mb33  MA  Wjjj  =  v^j 


<5> 


Define  w  =  <wi  wn  win> 

<6> 

and  v  ■  <vi  *n  vm> 

It  is  easy   to  verify  that  the  equations   <5>  are  equivalent 
to 

MA  W  MB  =  V  <7> 

Solution  of  <7>  can  be   accomplished  by  standard  Gaussian 

elimination  procedures.    Specifically,   the  following  steps 

are  required. 

1)  LDLT  factoring  of  MA  (n  x  n) . 

Forward  reduction  and  back-substitution  for  e  right- 
hand  side  vectors. 


ID 

(3)    LDLT  factoring ~ of  MB  ,(e  x,  e) 


Forward  reduction  and  back- substitution  for  n  right- 
hand  side  vectors. 

This  entire  process  is  economical  of  both  storage  and  arith- 
metic operations  because  of  the  tridiagonal  structure  of  MA 
and  MB. 


Boundary  Conditions 

The  cyclic  boundary  condition  is  implemented  by  repeating 
the  node  numbers  of  the  western  boundary  on  the  eastern 
boundary  as  shown  in  Fig.  1.  As  already  noted,  this 
accounts  for  nonzero  entries  in  the  upper  right-hand  and 
lower  left-hand  corners  of  MA. 

It  is  sometimes  required  to  impose  a  Dirichlet  boundary 
condition  on  the  southern  and  northern  boundaries  of  the 
region.  Specifically,  the  subvectors  w  and  w  of  the  solu- 
tion vector  w  are  prescribed.  To  implement  this  boundary 
condition  the  following  modifications  to  the  standard 
solution  procedure  are  required. 

In  the  n  x  e  matrix  V  on  the  right-hand  side  of  <7>  the 
first  and  last  columns  are  replaced  by  the  prescribed  bound- 
ary values  of  w,  i.e.,  put  v,  =  w,  and  ve  =  w  g.  Let 
X  =  W  MB  and  solve  the  system 

MA  X  =  V  <8> 

processing  successive  columns  of  V  in  standard  fashion,   but 
omitting  the  first   and  last  columns.    The   reduced  problem 


now  takes  the  form  W  MB  =  X  and  the  first  and  last  columns 
of  X  are  w,  and  we  ,  respectively.  Transposing  both  sides  of 
this  equation  gives 

MB  WT  =  XT  <9> 

where  WT  and  XT  are  the  respective  transposes  of  W  and  X 
(recall  that  MB  is  symmetric  and  is  thus  not  altered  on 
transposition).  Since  the  first  and  last  rows  of  WT  are 
known,  the  corresponding  scalar  equations  are  not  needed. 
Accordingly,  we  form  MB1  by  deleting  the  first  and  last  rows 
of  MB.  We  also  reduce  XT  to  XT1  by  omitting  the  first  and 
last  rows.   This  leaves  the  result 

MB1  WT  =  XT1  <10> 

or,  in  extenso,  this  takes  the  form  (for  n  =  3,  e  =  5) 


mb  2 1 

mb2  2 

mb  2  3 

0 

0 

mb  3  2 

mb  3  3 

mb  3  tt 

_  0 

0 

mbij  3 

mb  I*  k 

~k 

k 

k" 

0 

u 

u 

u 

"k   k   k~ 

0 

u 

u 

u 

= 

k   k   k 

mbi»  5_ 

u 

u 
k 

u 
k. 

_k   k   k_ 

(In  WT  and  XT1  the  elements  denoted  by  "k"  are  known  and 
those  denoted  by  "u"  are  unknown. )  This  equation  may  be  put 
in  standard  form  by  first  altering  the  first  row  of  XT1  by 
subtracting  mb21  times  the  corresponding  entries  in  the 
first  row  of  WT  and  altering  the  last  row  of  XT1  by  sub- 
tracting mbi,5  times  the  corresponding  entries  in  the  last 
row  of  WT.  Calling  the  new  right  hand  side  XT2  and  forming 
MB2  from  MB1  by  discarding  the  first  and  last  columns  and 
forming  WT1  by  discarding  the  first  and  last  rows  of  WT,  the 
result  is 

MB2  WT1  =  XT2  <11> 

Solution  of  <11>  is  carried  out  by  LDLT  factoring  of  MB2 , 
followed  by  forward  reduction  and  back  substitution. 


Floating  Point  Operations  and  Matrix  Storage  Requirements 

Presented  here  is  a  comparison  of  floating  point  opera- 
tion counts  and  matrix  storage  requirements  for  the  tensor 
product  scheme  and  three  widely-used  solution  algorithms  for 
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solving  <7>  or  its  equivalent  <1>.  Three  of  the  schemes 
take  advantage  of  symmetry  of  the  coefficient  matrices  and 
store  only  elements  on  or  above  the  principal  diagonal.  One 
of  these,  the  "band  solver"  (BAND),  places  these  elements  in 
a  rectangular  matrix  ne  x  r,  where  r  is  the  maximum  row 
length  of  the  upper  triangle  of  M.  The  "sky-line  solver" 
(SKY)  further  economizes  by  storing  only  that  part  of  the 
upper  triangle  beginning  at  the  diagonal  and  extending  to 
the  topmost  nonzero  element  of  each  column.  These  subvec- 
tors  are  assembled  into  a  single  vector.  This  scheme 
requires  an  additional  integer  address  vector  of  length 
ne  +  1.  The  remaining  algorithm,  "successive  over-relax- 
ation" (SOR),  is  iterative  rather  than  direct. 

In  most  applications  of  the  direct  solvers  the  number  of 
floating  point  operations  required  to  factor  the  coefficient 
matrix  into  LDLT  form  is  much  greater  than  those  required  to 
complete  the  process  of  finding  a  single  solution  vector  w 
corresponding  to  a  given  right-hand  side  vector  v  (forward 
reduction  and  back-substitution).  In  the  present  applica- 
tion, however,  the  latter  solution  process  must  be  carried 
out  17  times  for  each  time  step,  so  that  the  LDLT  factoring 
makes  a  negligible  contribution  to  the  total  computational 
expenditure.  Accordingly,  the  operations  required  for  fac- 
toring are  not  included  in  the  tabulation  below. 

In  the  following  table  the  results  given  for  the  number 
of  floating  point  operations  are  given  in  terms  of  the  grid 
parameters  n  and  e  (defined  in  Fig.  1).  One  multiplication 
(or  one  division)  plus  one  addition  (or  one  subtraction)  is 
counted  as  one  operation.  Exact  results  for  these  operation 
counts  would  take  the  form  of  a  polynomial  in  n  and  e.  Only 
the  highest  degree  terms  are  given  in  the  table.  Since  it 
is  not  possible  to  predict  the  number  of  iterations  per 
solution  when  using  SOR,  the  operation  count  given  for  that 
algorithm  is  for  a  single  iteration.  Also,  since  the  number 
of  storage  locations  required  for  SOR  coefficient  matrices 
is  highly  grid-dependent,  no  such  entry  is  given  for  SOR. 


TABLE  I.  Operation  Counts  and  Storage  Requirements 


Algorithm 

Number  of  Operations 

Number  of  Storage  Locations 

i 

per  Solution 

for  Coefficient  Matrices 

SOR 

10  en      (1) 

(2) 

SKY 

2  en2 

en2 

BAND 

4  en2 

2  en2 

TENSOR 

8  en 

3  n  +  4  e 

Notes:  1.  Number  of  operations  per  iteration. 

2.  Number  of  storage  locations  is  grid- dependent . 


Conclusion 


Close  comparison  of  operation  counts  and  storage  require- 
ments leads  to  the  conclusion  that  the  TENSOR  algorithm  is 
clearly  superior  to  the  SKY  and  BAND  algorithms.  The  com- 
parison with  SOR  is  not  as  clear-cut.  Considering,  however, 
that  the  operation  count  for  SOR  is  for  only  one  iteration, 
there  really  seems  to  be  little  doubt  that  TENSOR  is  the 
method  of  choice. 
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APPENDIX  A 
MATRICES  MA,  MB,  AND  THE  TENSOR  PRODUCT 


Symbols  a.  and  b.  which  appear  in  MA  and  MB  are  defined 
in  Fig.  1. 


MA   =  | 


(n    =    4) 


2(a<,+  aj  ai  0  a., 

ai        2(ai+  a2)  a2  0 

0  a2        2(a2+  a3)  a3 

a^  0  a3        2(a3+  a*) 


MB    = 


(e   =   3)      L 


2bi  bx  0 

bi    2(bx+  b2)       b2 

0  b2  2  b2 


The   tensor  product    of  matrices    C   and  D  may   be   represented 
in  block  partition   form  as 


C*D     = 


di  D     Ci2  D     ci3D 

c21 D     c22  D      c23  D 

_c31D      c32D      c33D 


where  the  c  .  are  the  elements  of  C.  Note  that,  if  C  and  D 
have  dimensions  r  x  s  and  t  x  u,  respectively,  the  tensor 
product  has  dimensions  rt  x  su. 
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APPENDIX  B 
FORTRAN  PROGRAM  LISTINGS 


FORTRAN  programs  for  the  implementation  of  TENSOR  are 
listed  here.  They  appear  in  the  form  of  subroutines  AMTRX2 , 
FACTOR,  BACKA,  and  BACKB  within  the  test  program  GAUSS3. 
(The  subroutines  FACTOR,  BACKA,  and  BACKB  are  adapted  from 
subroutine  COLSOL  of  Ref.  5.)  Also  included  are  G0G3 ,  an 
Exec  used  to  execute  GAUSS2  -  a  dimensioned  version  of 
GAUSS3,  and  CDIM,  an  Xedit  program  used  to  enter  the  dimen- 
sions . 
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Listing:   GAUSS 3  FORTRAN 

C   MAIN  PROGRAM   MASS  MATRIX  USING  TENSOR  PRODUCT  FACTORS 

C 

C    THIS  PROGRAM  IS  DESIGNED  TO  TEST  THE  SCHEME  (TENSOR) 

C    WHICH  RESOLVES  THE  MASS  MATRIX  INTO  A  TENSOR  PRODUCT  IN 

C    ORDER  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS  M  w  =  v.   THE 

C    SUBROUTINES  MAY  BE  INSERTED  IN  THE  PROGRAM  DEVISED  BY 

C   HINSMAN. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM 1 A/ NLAT, NLONG 

COMMON /CM8/A(Z1).B(Z1)  v    ,   % 

COMMON  AG(ZB) ,BG(ZC) ,GA(ZK) ,GB1(ZL) ,GB2(ZL) ,MA(ZM) , 
1MB(ZN)      ,   x 

DIMENSION  V(ZP) 

READ(5,  *) NLONG,  NLAT 

LATX=NLAT+1 

WRITEI6.1000) 
1000  , FORMAT (//,'    MASS  MATRIX  -  TENSOR  PRODUCT  RESOLUTION 

WRITE (6.1001) NLONG , NLAT 

REAd75.-)AaB 

WRITE? 6, 500) A 


WRITEi6,503)B 

FORMAT(/,'  B:  ' ,(24F3.0)) 
FORMAT (/  f  A:  \(24F3.0)) 
FORMAT T     NLONG  =  *I3, 


503 
500 

1001  FORMATC'    NLONG  =  *  ,  13  ,  '      NLAT  =',I3,/) 
C   CONSTRUCT  FACTORS,  GA,  GB1,  AND  GB2 ,  OF  MASS  MATRIX 

CALL  AMTRX2 
WRITE (6. 501) AG 

501  FORMAT?//         AG:       ' ,(12F4.1)) 
WRITE?6.504)BG 

504   FORMAtT;,'    BG:   ' .(12F4.1)) 
WRITEI6,i002)GA 

1002  FORMAtT/,'    GA' , / , (3X, 6F7 . 3 ) ) 
WRITE(6,i004)GBl 

1004  FORMAT^/.'    GB1'  ,  /  ,  (3X,  6F7  .  3  )  ) 
WRITE (6, 1005 )GB2 

1005  FORMAT?//  '    GB2 ' , / , (3X, 6F7 . 3) ) 
WRITE(6,l603)MA 

WRITE (6. 1006 )MB 

CU=GB1?3) 

CL=GB1U"LATX-1) 

K= (NLAT- 1) -NLONG 
C   IF  NDIR>0,  THERE  IS  A  DIRICHLET  BOUNDARY  CONDITION  ON 
C   NORTHERN  AND  SOUTHERN  BOUNDARIES. 

READ(5    *)NDIR,V 

WRITE?6\5025nDIR,V  . 

502  FORMAT?/, '  NDIR=  ',11,/,'  V:  f , 6F8 . 2 , / , (4X, 6F8 . 2) ) 
C   PERFORM  LDLT  FACTORING  OF  GA,  GB1,  AND  GB2 

CALL  FACTOR (GA, MA, NLONG) 

CALL  FACTOR (GB1, MB, LATX) 

CALL  FACTOR fGB2, MB, LATX) 

WRITE (6, 1002 )GA 

WRITE(6,1004)GB1 

WRITEI6,1005)GB2 
C   PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
C   FACTORS  OF  GA 

CALL  BACKA(GA,V,MA,NDIR) 
C   DIRICHLET  BOUNDARY  CONDITION  ON  NORTH  AND  SOUTH 
C   BOUNDARIES  ? 

IF(NDIR.GT.0)GO  TO  3 

WRITE?6,510W 
C   PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
C   FACTORS  OF  GB1 

CALL  BACKB(GB1,V,MB,NDIR) 

GO  TO  6 
C   CORRECT  RIGHT-HAND  SIDE  FOR  DIRICHLET  CONDITION 
3      DO  2  J=l, NLONG 

V 7 J+NLONG ) =V ( J+NLONG ) - CU*V ( J ) 
2     yiJ+K)=V?J+Kj-CL"V(J+NLAT"NLONG) 

WRITE (6, 5 10  )V 
C   PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
C   FACTORS  OF  GB2 

CALL  BACKB(GB2,V,MB,NDIR) 
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6  WRITE (6, 5 10 )V  , 

510   FORMAT (J  '  V:  ' , 6F8 . 2 , / , (4X, 6F8 . 2) ) 

C   READ  A  NEW  RIGHT-HAND  SIDE  AND  PERFORM  A  SECOND  SOLUTION 

READ(5  *)NDIR,V 

WRITE? 6\  502  )NDIR,V 

CALL  BACKA(GA,V,MA,NDIR) 

ifTndir.gt.ojgo  TO  7 

WRITE?6,510)V 

CALL  BACKB(GB1,V,MB,NDIR) 

GO  TO  16 

7  DO  5  J=l,NLONG 

V(J+NLONG)=V(J  +  NLONG)-CU*Y(.J) 

V(J+K)=vTj+Kj-CL-V(J+NLAT"NLONG) 

WRITE (6, 5 10 }V 

CALL  BACKB7gB2,V,MB,NDIR) 


5 


16 

1003 

1006 


WRITE 76. 5 10 )V 

FORMAT?/,     MA: ' ,2X,36I3) 

FORMAT(/  '    MB: '  2X  3613) 


STOP 

L  A  A  A  -•-  A  »•-  »•-  »•-  -'-  -•-  A  A  »»-  A  -.'-  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A 


^  **»  **  ^%"i%   Tji*  »i"  **  #**'%  V»  ^\1*   ******  *%   **  '***'%"**  **i"  **  *%  *i"  #*  W  *»  #*  **  *»  **  #%.  *i"  #%  *****  **  #%  #i"**  **  **  **%  **"****  **  *i*  #*  ** 

SUBROUTINE  FACTOR ( A, MAXA,NN) 
C 

C  .   -  -  INPUT  VARIABLES  -  - 

C  .  A(NWK}     =  STIFFNESS  MATRIX  STORED  IN  COMPACTED  FORM 

C  .  MAXA(NNP)  =  VECTOR  CONTAINING  ADDRESSES  OF  DIAGONAL 

C  .  ELEMENTS  OF  STIFFNESS  MATRIX  IN  A 

C  .  NN         =  NUMBER  OF  EQUATIONS 

C  .  NWK        =  NUMBER  OF  ELEMENTS  BELOW  SKYLINE 

C  .   -  -  OUTPUT  -  - 

C  .  A(NWK)     =  D  AND  L  -  FACTORS  OF  STIFFNESS  MATRIX 

C  .  .  . . 

IMPLICIT  REAL-8CA-H.O-Z) 

DIMENSION  A(1),MAXA(1) 

C      PERFORM  L*D*LT  FACTORIZATION  OF  STIFFNESS  MATRIX 


§ 


0     DO  140  N=1,NN 
KN=MAXA(N) 
KL=KN+1 

KU=MAXA(N+1)-1 
KH=KU-KL 
IF (KH) 110 ,90, 50 

IC  =  6 
KLT=KU 

DO  80  J=1,KH 
IC=IC+1 
KLT=KLT-1 
KI=MAXA(K) 
ND=MAXAfK+l)-KI-l 
IF(ND)80,80,60 
60     KK=MINO(IC,ND) 
C=0. 
DO  70  L=1,KK 

70   c=c+a(ki+l)*a(klt+l) 
a(klt)=a(klt)-c 

80     K=K+1 
90     K=N 

B=0. 

DO  100  KK=KL,KU 

K=K-1 


-J) 

100    A(KK)=C 

AIKN)=A(KN)-B 
110    IF ( A (KN 1)120.120, 140 
120    WRITE(l0UT,2d00)N,A(KN) 

STOP 
140    CONTINUE 

2000   FORMAT( / / » '  STOP  -  STIFFNESS  MATRIX  NOT  POSITIVE 

".ON 


1DEFINITET  ,//,'  NONPOSITIVE  PIVOT  FOR  EQUATK 
**&&•      M*OT  =  *,E20.12) 


END 
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A  «.•-  »•-  »*-  ■*•  *.'-  A  4*  ->~  -•-  •*•  4f  4m  -a.  -'-  *fr  *  *'-  -'-  *•*  -**  *fr  *  -**  **-  «'-  *  -*-  •"-  •*'*  *'«  -'-  ■**  -fr  -'-  -'-  -*'  -*~  *•'-  »'-  *  »'-  -*■  -*-  »'-  *'-  ••-  -*-  -*-  »'-  -*-  ••-  »'-  -'-  »'-  •*'-  »'-  »•- 


C*.*«  A  »»»  »*-  »'~  »»»  *.»-  »''-  »».»  V*  fcT*  *.'*  mfm  »»*  »**  •>■«  »'-  »*«•  »*-  »»-  *T*  »**  »**  -'-  »'-  *T*  ***  »**  *.**  »»~  «*'  »'-  «.*-  »'-  -.'*  «.»-  »**  -'-  «.'*  »*-  fc»*  *V  «A*  »'- 
«"t  •»  **  r\    *x  /»  *-%  **  <-»  **  #%  •>  «*'%««  ,r»  *-\     *x  *>  ^»  *»  *%    *^  »*  /*  »  *«  s>  s*.  '*  *>  *-i.  t^    r\    <-**»*»  *x  *■*    n*  -x  *\  /»  ,x 


-•..•. 


c  *■ 

C   THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 

C   SUBSTITUTION  USING  THE  FACTORS  OF  GA 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM 1A/ NLAT .NLONG 

DIMENSION  A(1),V(1),MAXA(1) 

C     REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 
C 

JMIN=1 

LATX=NLAT+1 

JMAX=LATX 
C   IS  THERE  A  DIRICHLET  BOUNDARY  CONDITION? 

IF7nDIR.LT. 1) GO  TO  140 
C   SKIP  NORTH  AND  SOUTH  BOUNDARIES 

JMIN=2 

JMAX=NLAT 
140   DO  240  J=JMIN.JMAX 
150    DO  180  N=l, NLONG 

KL=MAXA(N)+1 

KU=MAXA(N+1)-1 

IF(KU-KL)180,160,160 
160    K=N 

C=0. 

DO  170  KK=KL,KU 

K=K-  1 
170    C=C+A(KK)*V(K+(J-l)*NLONG) 

V{N+(J-1)*NLONG)=v(n+(J-1)*NLONG)-C 
180    CONTINUE 
C 

C     BACK- SUBSTITUTE 
C 

DO  200  N=l, NLONG 


K=MAXA(N) 
"111)* 


200    V(N+(J-l)*NLONG)=V(N+(J-l)*NLONG)/A(K) 
N=NL0NG 
DO  230  L=2, NLONG 
KL=MAXA(N)+1 

KU=MAXA(N+lI-l 

IF(KU-KL)230,210,210 
210    K=N 

DO  220  KK=KL,KU 

K=K-  1 
220    V(K+(J-l)'fNLONG)=V(K+(J-l)"NLONG)-A(KK)"V(N+(J-l)* 

1NLONG) 
230    N=N-1 
240    CONTINUE 

RETURN 

END 

c  .*.  ...  *     a,  ***i 

rj  ^\^%%%^\*\\c^\ *i\^%*}% ^>c^% *%'*%*}% ****** *i%**% '*£'*% ^%*}C^%'*\  ****** %**i\*yc**%^%  *+£*+%  *c*^*» "5% *c*iC  #c  vc*c  **  *%  *c  #»*»#%»»**#%#»*»'*  *% 

^*§HRQKHS1*§AC§1£4A 

C   THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 

C   SUBSTITUTION  USING  THE  FACTORS  OF  GB1  OR  GB2 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM 1A/ NLAT .NLONG 

DIMENSION  A(l) , V( 1) ,MAXA( 1) 

C      REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 
C 

LATX=NLAT+1 

NMIN=1 

NMAX = LATX 
C   IS  THERE  A  DIRICHLET  BOUNDARY  CONDITION? 

IFTnDIR.LT. l)GO  TO  50 
C   SKIP  NORTH  AND  SOUTH  BOUNDARIES 

NMIN=2 

NMAX = NLAT 
50     DO  240  J=l, NLONG 
150    DO  180  N=NMIN,NMAX 

KL=MAXA(N)+1 
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KU=MAXA(N+1)-1 

IF(KU-KL)180,160,160 
160   K=N 

C=0. 

DO  170  KK=KL,KU 

K=K-1 
170    C=C+A(KK)*V(J+(K-l)*NLONG) 

V(J+(N-lVA'NLONG)=v(j+(N-l)*NLONG)-C 
180    CONTINUE 
C 

C     BACK- SUBSTITUTE 
C 

DO  200  N=NMIN,NMAX 

K=MAXA(N) 

200   v  ( j  + (n-  1 ) *nlong ) =v ( j+ (n- 1 ) *nlong ) / a(k ) 
lmax=latx 

IF(NDIR.LT.l)GO  TO  205 

LMIN=3 

LMAX=NLAT 
205    N=LMAX 

DO  230  L=LMIN,LMAX 

KL=MAXA(N)+1 

KU=MAXA(N+lJ-l 

IF(KU-KL)230,210,210 
210   K=N 

DO  220  KK=KL,KU 

K=K-1 
220    V(J+(K-l)*NLONG)=V(J+(K-l)*NLONG)-A(KK)*V(J+(N-l)* 

1NLONG) 
230    N=N-1 
240    CONTINUE 

RETURN 

END 
C 

(~i       J^  .•..•.  .•. . .;.  .".  .'.  -'-  .'.  .•.  .'. .'.  .'.  J>.  -'..  -•-.  - '. . .'.  -'-  .'.  .•.  .'-  .'.  -•-  .' •  .'.  .'•  -'-  .'.  .'-  .'.  -'.  -'.  .•.  - '-  ■.•-  ."-  .'.  .'-  .*.  -•-  •.'-  ••»  .'.  .'»  .'.  .'.  .'-  •••  -'-  .'.  ••-  -••  -'-  -'-  -•-  -'. 

SUBROUTINE  AMTRX2 

C  THIS  SUBROUTINE  FORMS  THE  MASS  MATRIX  IN  THE  FORM  OF  A 

C  TENSOR  PRODUCT  OF  THE  GB  MATRIX  AND  THE  GA  MATRIX. 

C  THE  FIRST  OF  THESE  IS  NLAT  +  1  BY  NLAT  +  1, SYMMETRIC, 

C  AND  TRIDIAGONAL.   THE  SECOND  IS  NLONG  BY  NLONG, 

C  SYMMETRIC,  AND  TRIDIAGONAL  EXCEPT  FOR  SINGLE  ELEMENTS 

C  IN  UPPER  RIGHT  HAND  AND  LOWER  LEFT  HAND  CORNERS.  GB  IS 

C  STORED  IN  SKYLINE  VECTOR  FORM  (UPPER  TRIANGLE  WITH  SPACE 

C  FOR  FILL-IN)  AS  GB1  AND  GB2 .   THE  LATTER  VERSION  IMPOSES 

C  A  DIRICHLET  BOUNDARY  CONDITION  ON  THE  NORTH  AND  SOUTH 

C  BOUNDARIES.   GA  IS  ALSO  STORED  IN  SKYLINE  VECTOR  FORM. 

C  INTEGER  ADDRESS  VECTORS  MB  AND  MA  ARE  ALSO  GENERATED. 
C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON / CM 1A / NLAT , NLONG 


COMMON/ CM8/A(Z1].B(Z1) 
COMMON  AG(ZB),BG(Z<r 


1MB(ZN) 

:6n 


) ,GA(ZK) ,GB1(ZL) ,GB2(ZL) ,MA(ZM) , 


C   DIMENSION  BG (NLAT), AG (NLONG), GB1(2*NLAT- 1) , 
C    GA(3*NLONG-3],MA(NLONG+l),MB(NLAT+2) 

LATX=NLAT+1 
C    FIND  BG  =  (ELEMENT  HEIGHT) /6. 

DO  2  J=1,NLAT 
2      BG(J)=B(l+NLONG*(J-l))/3. 
C    GENERATE  GB1  AND  GB2 

GB1(1)=2.*BG(1) 

GB2(1)=1. 

DO  4  J=2,NLAT 

K=2*(J-l) 

GB1(K)=2.*(BG(J-1)+BG(J)) 

GBl(Ktl]=BG|j-l) 

GB 1 ( 2 -NLAT )  =  2 . *BG (NLAT ) 
GB1(2*NLAT+1)=BG(NLAT) 
GB2(3)=0. 
GB2(2"NLAT)=1. 
GB2(2*NLAT+1)=0. 
C    FIND  AG  =  (ELEMENT  WIDTH) /6. 
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GB2(K)=GB1(K 
4      GB2(K+l)=Gir 


DO  10  J=l,NLONG 
10    AG(Jl=A(J)/3. 
C    GENERATE  GA 

GA(l)=2.-(AG(l)+AG(NLONG)) 

DO  12  J=2,NLONG 

K=2*7J-1) 

GA7k)=2."(AG(J-1)+AG(J)) 
12     GA(K+1)=AG(J-1) 

Kl=2*NLONG 

K2=3*NLONG-4 

DO  14  K=K1,K2 
14     GA(K)=0. 

GA]3"NLONG-3)=AG(NLONG) 
C    GENERATE  DIRECTORY  VECTORS 

MB(1)=1 

D016  J=1.NLAT 
16     MB(J+1)=2*J 

MB(NLAT+2)=2*(NLAT+1) 

MA(1)=1 

DO  18  J=2,NLONG 
18    MA(Jj=2*p-l} 

MAiNLONG  + 1 )  =  3  * NLONG - 2 

RETURN 

END 
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Listing:   G0G3  EXEC 

ERASE  GAUSS 2  *  Al 

COPY  GAUSS 3  FORTRAN  Al  GAUSS2  = 

&STACK  CDIM 

&STACK  FILE 

X  GAUSS 2  FORTRAN 

FORTGI  GAUSS 2 

GLOBAL  TXTLIB  FORTMOD2  MOD2EEH 

FILEDEF  05  DISK 

LOAD  GAUSS 2  (START 


Listing:   CDIM  XEDIT 

SET  CMSTYPE  HT 
TOP 

*  ZB=NLONG 

C  /ZB/6/  *  * 
TOP 

*  ZC=NLAT 

C  /ZC/3/  *  * 
TOP 

*  Zl=NLONG*NLAT 
C  /Zl/18/  *  * 
TOP 

*  ZK=3*NLONG-3 
C  /ZK/15/  *  * 
TOP 
ZL=2*NLAT+1 

C  /ZL/7/  *  * 
TOP 

*  ZM=NLONG+l 

c  I  mm   *  * 

TOP 

*  ZN=NLAT+2 
C  /ZN/5/  *  * 
TOP 

*  ZP=NLONG*(NLAT+l) 
C  /ZP/24/  *V* 

TOP 

SET  CMSTYPE  RT 
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